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The problem of calculating real-time correlation functions is formulated in terms of an imaginary- 
time partial differential equation. The latter is solved analytically for the perturbed harmonic oscil- 
lator and compared with the known exact result. The first order approximation for the short-time 
propagator is derived and used for numerical solution of the equation by a Monte Carlo integration. 
In general, the method provides a reformulation of the dynamic sign problem, and is applicable to 
any two-time correlation function including single-particle, density-density, current- current, spin- 
spin, and others. The prospects of extending the technique onto multi-dimensional problems are 
discussed. 

PACS numbers: 02.70.Lq 



I. INTRODUCTION 

Quantum Monte Carlo (QMC) simulations proved very 
successful in studying many-body systems with realis- 
tic inter-particle potentials. Specific QMC methods have 
been developed for interacting electrons (for a recent re- 
view see jjj), liquid 4 He ||, alkali Bose-condensates in 
harmonic traps H , electron-proton plasma B and other 
systems. These methods yield accurate values of various 
physical properties, such as the cohesive energies of solids 
H or effective masses of defects || , which are in excellent 
agreement with experiment. In model problems includ- 
ing spin [Q, strong correlation ^j, and electron-phonon 
P,[lo| models, QMC usually handles much larger systems 
than other methods, treats equally well simple and com- 
plex interactions, and produces accurate results. 

A major difficulty faced by QMC is its inability to com- 
pute reliably dynamic properties of quantum mechanical 
systems. This is because QMC operates in imaginary 
time. For those quantities that do not directly trans- 
late in an imaginary-time language, the simulation noise 
quickly exceeds the signal rendering the calculation im- 
possible. This difficulty is known as the "dynamical sign 
problem" . It is usually dealt with by computing correla- 
tion functions in imaginary time first, and then continu- 
ing the results into the real-time domain. For the latter 
step, Pade approximants E ] , the maximum entropy [ pd| , 
singular- value decomposition and stochastic meth- 
ods are employed. Unfortunately, the analytic con- 
tinuation of noisy data is an ill-posed mathematical pro- 
cedure that does not lead to a unique solution. Some- 
times such an uncertainty results in spectacular failures 
of the reconstruction strategy. For instance, even the 
most advanced and popular method, the maximum en- 
tropy, cannot resolve the two peaks of the dynamic struc- 
ture factor in liquid helium [Q. This is a consequence of 
dealing with the mathematically ill-defined problem. An- 
other difficulty is assigning meaningful error bars to the 



continuation procedure. 

Thus the reconstruction procedure, while sometimes 
working reasonably well (especially in model systems), 
remain approximate, spoiling the rigor of QMC. What 
would improve the situation is the direct link between 
simulations and dynamic properties or rendering the 
problem mathematically well-posed. One such method 
was proposed by Mak and Egger who introduced "the 
multilevel blocking algorithm" to deal with the dynamic 
sign problem |l5f| . (A similar method was proposed for 
the fermion sign problem |l6|| .) Essentially, this was 
an attempt of straightforward evaluation of real-time 
rapidly oscillating path integrals that express the time 
dependent transition amplitudes. This method requires 
very large computer memory even for simple one-particle 
systems. Since the sign problem becomes more severe 
with increasing the size of the system, the blocking al- 
gorithm may not sustain the increase in the number of 
simulated particles. It leaves the stimulus for the search 
of other alternatives. 

In this paper, I propose another approach to alleviat- 
ing the dynamic sign problem. The calculation of real- 
time correlators is formulated as a mathematically well- 
defined initial value problem for the operator (ui — H) 2 
which is of the fourth order in spatial coordinates. The 
solution of the corresponding imaginary-time evolution 
equation yields directly the frequency-dependent spec- 
tral function, thereby eliminating the need for the time 
Fourier transformation, analytic continuation and other 
intermediate procedures. The basic idea is very general 
and is applicable to any quantum mechanical system and 
any two-time correlation function. It is described in Sec- 
tion ||. In Section III I apply the method to an exactly 
solvable case of disturbed one-dimensional harmonic os- 
cillator. Solving the new equation in more general cases is 
hard and require stochastic Monte Carlo methods. The 
arising difficulties and their possible solutions are dis- 
cussed in Sections |Tv| - |VIl| . 
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II. A DIFFERENTIAL EQUATION FOR THE 
SPECTRAL FUNCTION 

Consider a typical problem of many-body physics, that 
of finding the single-particle zero-temperature correlator 
in an interacting iV-particle system: 



tf(q,t) = (G| Cq (*) C t(0)|G) 

= e^ t (G| Cq e- m 4|G). 



(1) 



Here |G), Eg, and H are the ground state, its energy, and 
the full Hamiltonian of the system, respectively. The par- 
ticles could be either fermions or bosons. The operator c q 
creates an extra particle with momentum q, and operator 
c q destroys it. [It has to be emphasized that the single- 
particle momentum operators are chosen for definiteness 
only. In fact, everything that will be said and written 
about those operators will be equally valid for any pair 
of operators.] The notation Ti = 1 is used throughout the 
paper. K(q,t) describes the decay of \^f) = c q |G) with 
time due to 1^) not being an eigenstate of H. An equiv- 
alent amount of information is contained in the spectral 
density function A(q,u>) which is a Fourier transforma- 
tion of if (q, t): 



dt 



/OO 
K{q,t)e 
-OO 

= {2ir){G\c*5(u+E G -H) c q |G). 



(2) 



Compare the structures of Eqs. ([j]) and (g). The corre- 
lation function is complex and has an oscillating kernel 
e~ . The latter is the origin of the sign problem since 
big positive and negative contributions tend to cancel 
each other with a very small residue. On the contrary, 
A(q, to) is positive definite and has a non-oscillating ker- 
nel, a delta-function, which intuitively makes it a better 
computational object. There is another, physical reason 
to prefer the function A(q,ui) over K(q,t). The time- 
dependent properties themselves are rarely measured ex- 
perimentally. What is usually inferred from experiments 
are namely the energy-dependent spectral functions, or 
related to them response functions. Thus rather than 
compute Ks with subsequent numerical Fourier transfor- 
mations, it would be much better to be able to compute 
^4s directly. Besides, the main interest represent low- 
energy responses at small lus. This requires the knowl- 
edge of Ks at large times t that is precisely where the 
sign problem is worst. 

The above considerations suggest to compute A(q,oj) 
instead of K(q, t) jl?]]. However, while being formally 
correct, Eq. (Q) is impractical due to the presence of the 
operator function 8{uj + Eg — H). There is no way of 
replacing H with an explicit functional of the particle 
coordinates (which is required for numerical integration) 
other than using a different representation of the delta- 
function that would allow evaluation of the matrix el- 
ements of H in various powers. The common Fourier 



expansion of the delta-function would take the situation 
back to Eq. (||) and to the sign-problem. However, it is 
not necessary to use an integral representation. Instead, 
I use a limit representation of the delta-function, specifi- 
cally the limit of a gaussian. After that Eq. W) becomes 



A(%uj) = (2tt) lim (G|c q 



^ e - aiu+ E G -Hf c t | G) . (3) 



Taking the limit now amounts to solving the imaginary- 
time evolution equation 



dip 



Er 



h) i>, 



(4) 



supplemented with the initial condition ip(a = 0) = 

The equation ([1|) and representation (|3|) are the main 
results of the paper. They show that dynamic correla- 
tors can in principle be obtained by solving an imaginary- 
time partial differential equation. The equation is supple- 
mented by an initial condition and therefore it possesses 
a unique solution. The whole problem is thus mathemat- 
ically well-defined. The full computational procedure is 
as follows: (i) Obtain a ground state wave function |G) 
and energy Eq with whatever equilibrium method that 
works for the given system, (ii) Construct new state 
I*) = c q |G). (iii) Solve Eq. (§ using \^f) as an initial 
condition, (iv) Perform the a — > oo limit and compute 
the scalar product with (\&| to get the spectral density 
for given frequency uj. (v) Compute the response func- 
tions via fluctuation-dissipation and Kramers-Kronig re- 
lations. 

Let me discuss the general properties of Eq. (^) . First 
of all, this equation is of the first order in "time" vari- 
able a and fourth order in spatial coordinates because of 
the term H 2 . Mathematically, this is a parabolic equa- 
tion "correct in Petrovski sence" fl9[] . Its properties are 
quite different from the familiar second order parabolic 
(diffusion) equations. In particular, its Green function 
is no positive-definite but instead can have positive and 
negative regions. Because of that, the proposed method 
does not solve the sign problem. Rather, it transforms 
the dynamic, or "time" , sign problem into a "space" 
sign problem which resembles the fermion one. More 
on this in Section 



IV 



It is instructive to compare equa- 
tion (|4|) with the imaginary time Schrodinger equation: 
—dtp/dp = (H — Eg)iP- The latter's formal solution 



is ip(P) = e 



_ -p(H-E G 



V(0). It shows that in the j3 



limit, only the ground state survives because all the other 
states are exponentially suppressed. This well known fact 
is employed for instance in Projected Monte Carlo to ex- 
tract the ground state out of an arbitrary initial state 
-0(0). The crucial observation here is that the energies 
of other states are all greater than Eg which makes the 
exponent of the evolution operator be negative for all but 
the ground state. But what if one wants to retain not the 
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ground state but a state (or states) with energy (Eg+w)1 
In this case one would need a similar evolution operator 
with the exponent being a negative and even function 
of (Eg + u — H). The simplest function is quadratic 
one which results in precisely the exponential kernel of 
Eq. (|3|). Thus the physical meaning of the representation 
(|J) is that in the course of evolution only excited states 
with excitation energy 10 survive and all the other states 
with excitation energies above and below u> are exponen- 
tially suppressed. Notice also that while the variable (3 in 
the Schrodeinger equation is related to the inverse tem- 
perature, the variable a in Eq. (Q) has no clear physical 
meaning. The dimensionality of a is (energy) ~ 2 . Lastly, 
one may think that many other variations of Eqs. (^|) and 
(f|) are possible since there are infinitely many limit rep- 
resentation of the delta-function. This does not seem to 
be the case. All other common representation of 5(x) in- 
volve either non-analytical or non-polynomial functions 
(examples are e - "'^' and a / (1 + x 2 a 2 ) . Both cases make 
it difficult to interpret the resulting formulae in terms of 
a simple evolution equation. Comparable or even better 
representations may exist but I could not find any. 

I now rederive Eq. (^) differently, by demonstrating its 
equivalence to another definition of the spectral function. 
Let \(f> m ) be a set of eigenstates for the (N + l)-particle 
system, H\cj) m ) = E m \cj) m ). In general, |\E') is not an 
eigenstate of H but it is always expandable in \<j) m ): 



E 



\<Pm){<Pm\ci\G). 



(5) 



Upon substitution of the expansion in Eq. (|^) the Hamil- 
tonian in the exponent is replaced with E m . After that 
the a — > oo limit results in a delta- function 5(uj + Eq — 
E m ) and Eq. (|3j> becomes 

A(q,w) = (2Tr)J2(GK\<f> m )(<t>m\cl\G)S(cj + E G -E m ), 

m 

(6) 

which is the usual definition of A(q, ui) psfl . 

The formalism presented above is valid at zero temper- 
ature only. Its finite-temperature generalization is not 
straightforward and has to be addressed separately. The 
correlation function is defined as 



(7) 



where \n) is a complete set of the iV-particle system, 
inverse temperature, and Z = exp (—f3E n ) the parti- 
tion function. A new difficulty in comparison with Eq. (Q) 
is the second real-time operator e lHt which acts in a dif- 
ferent state space than e~ lHt . Therefore, a time Fourier 
transform of Eq. (j?j) does not lead to a convenient repre- 
sentation of the spectral density similar to Eq. (|^) . Put 
differently, e~ lHt cannot be replaced with a c-number be- 
cause it now acts not on a single state but on the whole 



set (n\ which is not supposed to be known in its entirety. 
Still, one can proceed further by formally separating the 
time variables in the two operators. Introducing an aux- 
iliary delta-function 1 = J dt'S(t — t'), replacing t — ► t' 
in e~ lHt , using 5(t — t') = J ||e _le ^ _t \ and integrating 
over t' one obtains 

K(q,t) = 
1 

- dee- iet ^{n\e lHt Cci 5{e- H)c^e- pB \n). (8) 

~ °° n 

A Fourier transformation yields for the the spectral den- 
sity 

n poo 

— / deJ2(n\ 5(lu-e + H) c q 5(e - H) c\e- 0H \n) 

J — oo m 



27r fai [et2 

— lim lim w — a — 

Z ai^oo 02^00 y 7f V 7T 



de 



X (n| e -°*i«-e+Hf Cq e - a2 (e-Hf c y,-l3H |^ 



(9) 



There are three imaginary-time evolution processes in the 
last equation. First one, starting with an arbitrary state 
|n), during time (3, and under operator H. Second one, 
starting from the finite state of the first process plus an 
extra particle with momentum q, during infinite time, 
and under operator (e — H) 2 . Third one, starting with 
the final state of the second process minus a particle with 
momentum q, during infinite time, and under operator 
(oj — e + H) 2 . In the end, the scalar product with (n\ and 
integration over the real variable e have to be done. In 
the zero-temperature limit, the operator e~@ H projects 
out all but the ground states, e~^ H \n) = e~ l3EG \G)5 n G, 
and then the factor e~^ E ° cancels Z. After that, the 
leftmost H in Eq. (|^) is replaced with Eg and the limit 
in ai and integration over e are easily performed, thereby 
reducing Eq. (|^) to the previous expression ([}]). 

The finite-temperature generalization is achieved at 
the expense of the second auxiliary evolution and an ex- 
tra integration. This makes any practical realization of 
the method much more difficult than at zero tempera- 
ture. Nevertheless, Eq. (|^) is exact and contains only 
real exponentials. 

In conclusion of this section, let me reiterate that 
nowhere the specifics of the single particle operators 
(fermionic or bosonic) have been used. Everything that 
has been said about c q and c q is equally applicable to 
any pair of operators 0\ and Oi- A generic expression 
for the zero-temperature spectral density is a straightfor- 
ward generalization of Eq. (0) : 



A 6 1 6 2 (") = ( 2 



tt) lim (G\6 2 \[^-e- a ^ +EG - H ^ 6 X \G). 



(10) 
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There are no restrictions on the choice of the operators 
Oi and O2. They could be single-particle, two-particle, 
spin, or other. In the case of the density or current oper- 



ators, Eq. (10) opens the exciting possibility to compute 
dynamic structure factors and conductivities from first 
principles. 



III. PERTURBED HARMONIC OSCILLATOR: 
ANALYTICAL SOLUTION 

The purpose of this Section is to solve a non-trivial ex- 
ample to demonstrate Eqs. (|J) and (|J) at work. Consider 
a harmonic oscillator with frequency Q and mass M. It 
can be either free of subjected to a constant force of mag- 
nitude gV2Alhtt 3 , depending on whether an additional 
particle is present in the vicinity of the oscillator or not. 
The Hamiltonian reads 



H , 



1 d 2 
'2'dx 2 



f 



a/2. 



,gx c 1 c 



(11) 
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FIG. 1. The wave function ©-([[I) at different "times" 
a. g — 1.5 and uj = 2.0 £1. Note how negative pockets develop 
with time. After a = 5.0 the solution changes only slightly. 



where the coordinate x is measured in units of \Jfij (Af£l) 
and g is the dimensionless coupling constant. This 
Hamiltonian is known from the independent boson model 
p8[ . It can also be thought of as the on-site Holstein 
polaron pc| ]. Let the oscillator be in its ground state 
and the particle absent at t < 0. The state of the sys- 
tem is \G) = lOpartjOz) = |0 part , -^ji^ x ' 1 ^) with energy 
Eq = ifi. At time the particle is placed on the oscilla- 
tor creating state \i/)(x, a = 0)) = |l par t, X ). At this time 
the oscillator begins to experience the force. At time t 
the particle is removed and the oscillator is left free but 
in one of the excited states \m x )- So between and t the 
Hamiltonian is 



V2. 



gx 



(12) 



Our objective is to compute the spectral function A(u) 
of operators cf and c. According to the general scheme, 
one requires to solve equation (^) with H x from (|l2] ) 
and initial condition |0 X ). This can be done by intro- 
ducing first a new variable y — x — \/2<? and then bosonic 
operators b y and b y . The transformed Hamiltonian is 
H y = Q(b y b y — g 2 + \). After expansion of the solu- 



tion in the eigenstates of b y b y , \tp) = Yll 
Eq. (0) becomes: 



l {a)\m y 



m=0 



da m (a) 



da 



m,y) =y^ a m (a)(uj + g 2 fl - £lb\by) 2 \m y ) 

m=0 

00 

a m (a)(uj + g 2 fl — ilm) 2 \m y ). 



(13) 



m=0 



Equations with different m are separated and easily in- 
tegrated yielding 



m=0 



1 



7T l/iy/2 m m \ 



H m (x~V2g)e * 



(14) 



where H m (x) are the Hermite polynomials. The co- 
efficients a m (0) are found from the initial condition 
ip(x,0) =7T- 1 / 4 exp (-x 2 /2): 



Om(0) = (m,l0 a ) = ( i^L e -- 



(15) 



V2" I m! 

Figure |l| shows the solution (14)-(fj"5|) as a function 
of "time" a. Note the appearance of negative regions. 
Finally, the scalar product (0 x \ip(x, a)) yields the "en- 
dependent spectral density" 



A(u,a) 



(2ir)e- 9 ' 



m=0 



-a(ui+g 2 Q-nm) 2 



(16) 



Figure |] illustrates the changes in A(u>, a) as a progresses 
from to 00. The solution evolves from the identical zero 
through a smooth single peak to a highly non-analytical 
grid of delta- functions. The a — ► 00 limit of Eq. ( |l6| ) 
results in 



A{u) = (2ir)e~ 9 ' 



00 2m 

^ ^5{u) + g 2 Q. - Om), 

rn—0 



m 



(17) 



which is known to be the correct spectral density |lq ] 
(chapter 4). 

Two comments are appropriate at this point. The first 
one concerns the solvability of the equation (Q). In the 
above example, the equation was solvable because the os- 
cillator before, during, and after the perturbation was an 
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G(X',X;a,e) = (X.'\e- a(e - H)2 |X) 



(18) 




2.0 4.0 6.0 8.0 

FIG. 2. The spectral function @ of the perturbed 
harmonic oscillator for different values of the parameter a. 
g = 1.5. The first (leftmost) peak appears at u) = —g 2 £l. The 
area under all the curves is 2ix. This is how the correct result 
would emerge in the course of numerical solution of Eqs. ([]) 
and (g). 



integrable system which eigenfunctions were well known. 
Surely, for any such a system where all the eigenfunctions 
are known, the equation (jj) would be solvable as easily as 
for the oscillator, and the results would be expressible in 
the form of a spectral expansion similar to Eqs. jl^)-(16). 
The second observation is about the differences between 
systems with discrete and continuous spectrum. In the 
case of discrete spectrum, as in the above example, the 
spectral function must have the form of a grid of delta- 
functions, see Eq. (^]). This is realized in the present 
formalism in the following way. For certain frequencies 
u>, the solution of equation (Ji|) has a finite overlap with 
the final state (G|c q even in the limit a — ► oo. Then the 
factor ^/a results in infinite limit values of the spectral 
function for these resonant frequencies. For systems with 
continuous spectra, including all real systems, the overlap 
and consequently the solution of (^) must decay as 1 / y/ot 
in the a — > oo limit, so that the product with the factor 
yfa remains constant. If such an asymptotic regime is 
observed at some large a, then this is an indication that 
the solution process can be stopped. 



IV. THE SHORT-TIME PROPAGATOR 

Equation (||) can be solved analytically in a hand- 
ful of cases. In many-body problems, the equation be- 
comes multi-dimensional and therefore requires stochas- 
tic methods of solution. The central task here is the 
calculation of the Green's function (GF) of the evolution 
operator: 



where X = {xi} and X' = {x'j} are two sets of the sys- 
tem's coordinates. The number of coordinates is n = Nd, 
N being the number of particles and d the dimensionality 
of space. I have also used the notation e to represent the 
sum uj + Eg- Let us find GF for the system of free dis- 
tinguishable particles with unit mass which is described 
by Hamiltonian 



H = T 



2 ^ ftr? 

l—l L 



-V 

2 



(19) 



where V„ is the n-dimensional gradient operator. Intro- 
ducing the momentum states |K) via 



dK 



iKX 



K) 



(20) 



and using (K'|K) = (27r)"(5(K — K') one obtains from 



G (X',X;a,£) = 



dK 

(27T)"' 



iK(X-X') e -a( E -J§-) 2 _ ^1) 



The last expression satisfy the initial condition 
Go(X',X;a = 0) = 6(X - X'). As expected from the 
translational invariance, GF is a function of the differ- 
ence (X — X'). The rotational invariance of space also 
suggests that GF must be the function of the modulus 
|X — X'| = R only. This is indeed the case which allows 
the reduction of Eq. ( |2l| ) to a one-dimensional integral 
(technical details are not important for the purposes of 
this paper): 



G {R;a,s) 



2tt- 

(2tt)" 



dkk ArX i e ~ a(£ ~ 



(22) 



where J v {x) is the Bessel function. In one dimension, 
n = 1, it follows from Eq. (|22|) or directly from Eq. (|2ll ) 
that 



dk 
2^ 



COs[fc(lE — x')] i 



(23) 



The last equation is a convenient place to compare the 
properties of the imaginary time Schrodinger equation 
(or diffusion equation) with those of Eq. (^) . In the for- 
mer case, the corresponding expression for Gq =1 would 
have just k in place of (e — -^-). Then, the integral would 
easily be performed leading to the positive-definite GF 
of the diffusion equation. In the present case, despite 
the integral cannot be performed analytically, it is easy 
to see that GF is not positive definite. Indeed, take the 
limit of large a. The intergal is dominated by the points 
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FIG. 3. The free one-dimensional Gren's function (|23|) 
for two values of the time interval a and three values of the 
excitation energy e. The oscillations grow with both a and e. 

k = ±^/2e and therefore as a function of la; — x'\ oscillates 



with period (27r)/v2e. This is the source of the negative 
regions and of the "space" sign problem. Notice that the 
period is large for small to and small for large u>. Thus 
the sign problem is expected to be less severe in the phys- 
ically interesting regions of low-energy excitations. This 
trend is encouraging and it is opposite to what happens 
in the time formulation, as discussed in Section [n]. For 
the purposes of this paper it is more interesting to study 
the case of small a and large distances \x— x'\. A rigorous 
analysis is quite involved, let me only mention that the 
amplitude of the GF decays exponentially in this limit: 



\G% =1 (x', x;a)\< C x exp (-C 2 \x - x'\ 4 / 3 ) , 



(24) 



3| shows the results of 
for several a and e. 



where G X ,G 2 > [19|. Figure 
numerical integration of Eq. (E3 
Notice that the number of oscillations increases with the 
excitation energy, a seems to affect the amplitude of 
the oscillations rather than their period. In conclusion 
of this part one should mention that the existence of the 
negative regions in the Green's function was shown to be 
a general property of differential operators with the order 
of spatial derivatives higher than 2 |H[ . 

Let me turn to the more difficult and more interesting 
case of interacting quantum mechanical system. The full 
Hamiltonian now consists of the kinetic and potential 
parts, H = T + V . The potential energy V may include 
single-body as well as two-body terms. As it is always 
done in Monte Carlo studies, the long-time propagator 
( |l8| ) can be expressed as a multiple integral of the short- 
time propagator 



G(X',X;a) = 
J dXi . . . dX L G(X', X L ; Aa)... G(X l5 X; Aa). (25) 

The number of additional time slices L must be cho- 
sen large so that the new time interval A a = ^"b; 
gets small enough to allow accurate approximations of 
G(Xj+i, Xj| Aa). The main difficulty in constructing 
those approximations is the non-commutativity of the 
kinetic and potential operators T and V. This difficulty 
is usually overcome by the Trotter-Suzuki decomposition 
of exponential operators, see, e.g., [E2j. However, in the 
present case one is faced with a completely new situation. 
The exponent of the evolution operator [the right hand 
side of Eq. (0)] contains not only sums but also products 
of non-commuting operators, namely TV and VT. Thus 
all the recipies have to be devised anew. In the first order 
in Acv it can be done in the following way. Expand the 
exponential in Eq. nUty: 

G(X l+1 ,X;Aa,e) = (26) 
(X, i+1 | {l - Aa(e -f- V)(e -f-V) + o(Aa)} |X). 

In the term linear in Aa, the potential operator in the 
first parantheses can be replaced with its value at the 
final configuration: V —* V^X,-_|_i) = V^+i. Similarly, 
the potential operator in the second parantheses can be 
replaced with its value at the initial configuration: V — > 
V(Xj) = Vi. Exponentiating back, one obtains 

G (Xi+i,Xi; Aa,e) 

W <X +1 \ e -^(e-V i+1 -f )(e- Vi -f )| Xi) 



= pttW+i-v*: 



i+i|e 



-Aa E -i(V i+ i+Vi)-f 



(27) 



e 4 



■ (v i+1 -v{y 



G ( X i+ i,Xj; Aa,e - gO^+i 



Vi) 



The last line expresses the short-time propagator of the 
interacting system via the short-time propagator of the 
non-interacting system Go taken at a shifted energy ar- 
gument. This result together with Eq. ( p2| ) lead to the 
following remarkable conclusion: even in a many-body 
interacting system calculation of the short-time propaga- 
tor requires knowledge of the potential energy at the end 
points and only one one- dimensional integration. More- 
over, since the integral in Eq. ( p2] ) is a universal function 
for given n and a, it can be easily tabulated for all re- 
quired values of s and |X^ + i — Xj|. 

Of fundamental importance is the question whether 
the higher order approximations could be derived. Pre- 
liminary studies suggest a positive answer albeit with 
much more cumbersome expressions than Eq. • Con- 
tinuing the expansion in Eq. ( p6| ) further one finds in the 
next term already four factors (e — T — V). Only two 
operators V can be replaced with numbers immediately, 
the other two have to be commuted through Ts. Such 



G 



a process generates multiple additional terms that con- 
tain (V n V)s and (V^V)s. I do not show specific details 
here but it is clear that the second order approximation 
would be much harder to evaluate numerically than the 
first order one, especially in the multi-dimensional situa- 
tion. Obviously, a special study is needed to investigate 
whether such a complexity is worth the gained accuracy 
and the consequent increase of the time step Aa. 



V. INTEGRATION 

With the short-time propagator known for given ex- 
citation energy u> and time step Aa, the "a-dependent 
spectral function" is given by the multiple integral 

A(a) = ( 27r )y^ J rfXdXi . . . dX L dX' x 

$(X') G(X', XiiAo)... G(Xi, X; Aa) tt(X), (28) 

where VP and $ are the initial and final states of the 
system. (Recall that ^> and $ may not be the same, 
this depends upon the choice of the operators 0\ and O2 
studied.) While trying to evaluate the above integral one 
faces an immediate difficulty of infinite limits of integra- 
tion. It can be resolved by introducing auxiliary normal- 
ized probability densities -P(X) with properties P(X) > 
and J dXP(X) = 1. Then the last equation is rewritten 
identically as follows 



A(a) 
2tt 



J dXdX 1 ...dX'W P'(X') ... Pi (Xi)P(X) 



TT JdXdX 1 ...dX'P>(X>).. 



.Pl(X!)P(X) 

(29) 



w 



_ $(X')G(X', X L ; Aa)... G(X l5 X; Aa) *(X) 



P'(X')P L (X L )...P 1 (X 1 )P(X) 



(30) 



As a result of this transformation, the product of func- 
tions P can serve as the weight function. The explicit 
form of functions P is completely at our disposal. It 
should be chosen to minimize the probability of large 
deviations of Xs. The functions P need not be all the 
same. Each P is absolutely arbitrary and independent 
of the other as long as the main conditions of positiv- 
ity and normalization are fulfilled. In short, functions 



P serve to optimize stochastic integration of Eq. (29). 
Note also that the absolute values of the Green's func- 
tions |G| cannot be used as auxiliary probability densities 
because their normalization is not known a priori. [One 
may think that one could get away from this obstacle 
by simply assuming some normalization factor and then 
restoring it afterwards from the overall normalization of 
the spectral function. While this can work in principle, 
this would require a very accurate determination of A 
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FIG. 4. The Green function @ for a = 2.0 fi~ 2 and 
uj = -2.75 n (top left), w = —1.75 Q (top right), uj = -0.25 CI 
(bottom left), and uj = 1.25 (bottom right). The number 
of nodes increases with excitation energy. The two horizontal 
axes are x and x' . 



even at those energies one is not interested in. Besides, 
not every spectral function satisfies a normalization sum 
rule.] 

The actual measured quantity is the ratio W from 
Eq. (|30|). Since Gs have negative regions their product 
can be both positive and negative. The sign alternation 
of W is going to be the main technical difficulty in im- 
plementing the present idea in practice. 



VI. PERTURBED HARMONIC OSCILLATOR: 
NUMERICAL SOLUTION 

The results obtained in the last two sections will now 
be applied to the harmonic oscillator example discussed 
above. Two issues are going to be addressed specifically: 

(i) the accuracy of the first-order approximation ( p7| ) and 

(ii) stability of the integration in (p9|)-(|30|). 

The exact oscillator propagator for any time interval 
a is given by 



G(x , x; a, uj) 



1 



.i(x'-V2g)%-i(x-V2 ff ) 2 



(31) 



OO 

2 m m! 

m— 



H m { x' -V2g)H m (x - V2g)t 



-a(u>+g Cl-tlm) 2 



This function is shown in Figure I for a = 2.0 0~ 2 and 
several excitation energies uj. Notice how the number of 
sign changes increases with u>. An analogy seems to exist 
between this fact and the increase in the number of nodes 
of the eigenstates of the time-independent Schrodinger 
equation. 

One can use the exact solution ( |3l| ) to check the ac- 
curacy of the approximate representation (E7T). Due to 



7 



the one-dimensional character of the problem the check 
can be very effectively done as follows. The approxi- 
mate short-time propagator is computed for some small 
time step Act. Then the propagator corresponding to 
the double time G(2Act) is found by convolving two 
G(Aa)s. If G(x' , x; ct)s are represented by matrices with 
matrix elements evaluated on a uniform mesh (x^Xj) 
then the convolution is equivalent to squaring of matrix 
G(x' i ,x j ;Aa) = Gy(Aa): 

Gy(2Aa) = Ax G 2 ,- (Act), (32) 

where Ax is the step of the coordinate mesh. Repeat- 
ing the squaring n times one arrives at the propagator 
corresponding to the large time interval 2™Act, which is 
compared with the exact result. In this manner it was 
established that the rule ( p7| ) is not as accurate as the 
one without the exponential factor exp[—p(Vj+i — Vi) 2 ]- 
For instance, at Act = 0.05 5!~ 2 the expression with- 
out the factor remains within 10% from the exact so- 
lution even after 9 squaring which corresponds to 512 
time steps or total time a = 25.6 i7~ 2 , in a wide range 
of energies —2.75 < u)/Ct < 2.75. At the same time, 
the formula with the exponential factor demonstrated 
deviations ~ 100% already at 7 squaring, i.e., at times 
ct w 6.4 Q~ 2 . At larger Act sa 0.2-0.4 fi" 2 , both versions 
were bad, Eq. (|27]) overestimating and the one without 
the factor underestimating the correct values of the long- 
time propagator. The best results were obtained with 
the factor but with the coefficient in the exponent being 
w 0.10 instead of 0.25 as in Eq. @. 

It is not the purpose of the paper to perform a detailed 
examination of this issue. Given the complex character 
of the new evolution operator, a complete analysis could 
be quite involved. A possible reason for the described 
inconsistency may be a bigger than expected role of the 
higher-order terms neglected in deriving approximation 
(p7j). A better understanding of this role requires ex- 
plicit derivation of the second and possibly third-order 
approximations and their thorough comparison. This is 
a difficult task that warrants a separate study. 

Let me now present the results of numerical integra- 
tion of Eqs. (^9"1)-(|30|). I have used gaussian functions 
P{x) and the simplest "brute force" method of integra- 
tion when uncorrelated sequence of configurations {x} is 
generated and the quantity W measured and accumu- 
lated in a straightforward manner. A crucial parameter 
in the process is the number of time intervals or the total 
number of integrations. As long as the number of inte- 
grations is less than about 20 the sign problem does not 
manifest itself and the simulations produce accurate re- 
sults with controllable error bars. An example is shown 
in Figure || where the exact propagator is used. The 
multipeak structure of the spectral function is clearly re- 
solved. Similar encouraging results have been obtained 
for other time steps and number of steps, for instance 



5.0 




-4.0 -2.0 0.0 2.0 4.0 

FIG. 5. Results of Monte Carlo integration of 

Eqs. (^)-(^) for the disturbed harmonic oscillator. The time 
step Act = 2.0 fi -2 and the exact propagator ( |3l| ) are used, 
see also Figure n. The full time interval is a = 20.0 fT 2 . The 
total number of integrations is 11. The number of measure- 
ments is 10 6 . 

Act = 0.8£T 2 ,iV = 16; Act = 1.0rr 2 ,N = 20; and so 
on. However, when the number of integrations exceeds 
20, the sign problem sets in and the simulations cannot 
be completed. 

The results of Monte Carlo integration with the ap- 
proximate propagator ( p7j ) are shown in Figure §. In this 
example, the numerical coefficient in the exponent of the 
exponential factor is taken to be 0.10. The number of 
time steps is 16 and the total "time" is ct = 6.4 fl~ 2 . 
Error bars have been estimated from the results of 5 
successive runs of 10 7 measurements. The agreement 
between the Monte Carlo data and the exact spectral 
function is good, given the simplicity of the numerical 
method used. In this respect several comments are ap- 
propriate. First of all, the perturbed harmonic oscillator 
is not a simple system as far as its spectral function is 
concerned. The latter has a multi-peak structure and is a 
hard task for any method. It is known that the inability 
to effectively resolve multiple peaks is the main short- 
coming of the mainstream reconstruction methods. The 
very fact that the present method is able to resolve mul- 
tiple peaks is very encouraging. This is a consequence of 
being a direct method: there is no difference between the 
smooth and the peaked function as long as the value for 
the given energy point is computed independently and on 
the same footing with other points. Secondly, it has not 
been the purpose of this paper to provide a full-scale nu- 
merical analysis of this example. The main interest lies 
in multi-dimensional interacting systems which neverthe- 
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VII. CONCLUSIONS 




-4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 



FIG. 6. Results of Monte Carlo integration of 

Eqs. (^)-(^o|) for the disturbed harmonic oscillator. The time 
step Aa = 0.4 Q.~ 2 and the approximate propagator ( p7j ) are 
used. The full time interval is a = 6.4 SI -2 . The total number 
of integrations is 17. (This is the reason for larger error bars 
as compared with Figure The number of measurements is 
10 7 . 



less may have simpler structure functions. There is not 
much sence in perfecting the solution when there is no 
guarantee that the tricks developed will work in the phys- 
ically more interesting cases. The results given above 
serve primarily as an illustration to the method and a 
demonstration of its viability. Thirdly, and perhaps most 
importantly, the "brute force" method of integration is 
not the right way to solve the equation (^). The incredi- 
ble progress of Monte Carlo methods in the past decades 
has been largely due to the "smart" ways of performing 
multiple integrations such as the Metropolis algorithm 
in statistical applications or the importance sampling in 
fermionic applications. In relation to the present paper, 
the recent exact fermionic algorithm of Kalos and Ped- 
eriva looks particularly interesting. They developed 
a Diffusion Monte Carlo algorithm for finding the ground 
state of the many-body fermionic Schrodinger equation 
with no approximation. Their method can handle wave 
functions with positive and negative regions, so it effec- 
tively solves the fermionic sign problem. Precisely an 
algorithm of this kind is required for solving equation 
(f|). It is too early to judge whether these ideas could be 
directly applied here but this is definitely one of possible 
directions to search for an effective method of integration. 



This paper describes a novel approach to direct 
stochastic calculation of dynamic correlation functions of 
quantum mechanical systems. It was inspired by the dif- 
ficulties of the traditional reconstruction methods orig- 
inating from the absence of a mathematically rigorous 
formulation of the analytical continuation procedure. I 
believe that at present the quest for an effective way 
of computing the dynamic correlators directly in real 
time/frequency domain is fully justified, and it should 
continue until any of the two approaches results in an 
effective and satisfactory numerical algorithm. 

Let me summarize how the present paper contributes 
to this process. It has been proposed to make the spectral 
function A(u>) rather than the time-dependent correlator 
Kit), the main object of the computational effort. This 
is advantageous because it directly leads to experimen- 
tally accessible properties and avoids the need for addi- 
tional procedures such as Fourier transformations. The 
calculation of the spectral function is made possible by 
replacing the operator delta- function with the limit of a 
gaussian and by interpreting the latter as a long-time so- 
lution of the evolution equation (||). Thus the process 
of finding A(u>) is formulated as an initial value problem 
for a fourth-order differential operator. Such a problem 
is mathematically well-defined and has a unique solution. 
One "only" requires to solve the equation in the limit of 
infinite time a. The solutions can be classified with re- 
spect to their limit behavior. Some will remain constant 
which will lead to infinite values of the spectral function 
(apparently at discrete energies only). Other solutions 
should decay as oc ^= and lead to finite A{w). Third class 
comprises solutions decaying faster than that (probably 
exponentially) and resulting in vanishing A(lu). In the 
end, taking a scalar product with a function of interest 
is required. The whole sequence should be repeated for 
all relevant excitation energies. The new formulation has 
been shown to be equivalent to the standard theory and 
definitions of the spectral function. It has also been gen- 
eralized to finite temperatures. Another positive feature 
of the method is its universality. It imposes absolutely no 
restrictions on the nature of the operators studied. They 
could be bosonic or fermionic, one-particle, two-particle, 
or even more complex, current, density, spin, and so on. 
Also the spectral functions are not required to possess 
any special properties such as non-negativity, additional 
symmetries, and sum rules. In short, the theory is sim- 
ple, mathematically well-defined, and universal. 

The practical usefulness of the proposed method will 
be very much dependent on the availability of stochas- 
tic methods that could solve equation (^) in multi- 
dimensional cases. Since the propagator necessarily has 
negative regions, handling of sign-alternating functions is 



required. An example considered in Section VI demon 
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strates that this is achievable. In broader terms, the 
situation bears strong similarities with Diffusion Monte 
Carlo studies of systems of interacting fermions, where 
a significant progress was reported p3 24 1. Additional 
work is needed to develop similar approaches for equa- 
tion (Q). If successful, the effort will bring direct ways of 
first-principle Monte Carlo studies of dynamic character- 
istics of quantum mechanical systems. 
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